{smcl}
{com}{sf}{ul off}{txt}{.-}
      name:  {res}<unnamed>
       {txt}log:  {res}C:\Users\Carolina\Downloads\reuse the seed\Tables 1-4 log_final.smcl
  {txt}log type:  {res}smcl
 {txt}opened on:  {res}29 Jun 2021, 11:16:26

{com}. do "C:\Users\Carolina\AppData\Local\Temp\STD00000000.tmp"
{txt}
{com}. ***************************************************************************
. ***************************************************************************
. **Replication files for Enns, Moehlecke and Wlezien: "Detecting True 
. **Relationships in Time Series Data with Different Orders of Integration"**
. **Tables 1-4 on main paper*************************************************
. **Last run: Jun 27th, 2021************************************************** 
. ***************************************************************************
. ***************************************************************************
. 
. version 14.2
{txt}
{com}. 
. *Note: we use the same seed throughout the simulations so that results are 
. *directly comparable. We reset the seed for each simulation.
. 
. 
. *seed selected based on random.org, b/t 1 and 1000000000
. set seed 457090451
{txt}
{com}. 
. 
. *************************************************************************
. *************************************************************************
. ****************************TABLE 1**************************************
. *************************************************************************
. *************************************************************************
. 
. 
. **********************************
. *Table 1 layout
. **********************************
. *generate excel file to store Table 1 output
. putexcel set table1.xlsx, sheet(sheet1) replace
{res}{txt}Note: file will be replaced when the first {cmd:putexcel} command is issued

{com}. *add labels to table
. *row labels
. putexcel A4 = "$\hat{c -(}\alpha{c )-}_1$"
{res}{txt}file table1.xlsx saved

{com}. putexcel A5 = "$\hat{c -(}\alpha^*{c )-}_1$"
{res}{txt}file table1.xlsx saved

{com}. putexcel A6 = "$\hat{c -(}\beta{c )-}_1$"
{res}{txt}file table1.xlsx saved

{com}. putexcel A7 = "$\hat{c -(}\beta{c )-}_2$"
{res}{txt}file table1.xlsx saved

{com}. putexcel A8 = "$\hat{c -(}\beta^*{c )-}_2$"
{res}{txt}file table1.xlsx saved

{com}. 
. *column labels
. putexcel B1 = "ADL"
{res}{txt}file table1.xlsx saved

{com}. putexcel B2 = "$\rho = 0.2$"
{res}{txt}file table1.xlsx saved

{com}. putexcel B3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel C3 = "%"
{res}{txt}file table1.xlsx saved

{com}. putexcel D2 = "$\rho = 0.5$"
{res}{txt}file table1.xlsx saved

{com}. putexcel D3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel E3 = "%"
{res}{txt}file table1.xlsx saved

{com}. putexcel F2 = "$\rho = 0.8$"
{res}{txt}file table1.xlsx saved

{com}. putexcel F3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel G3 = "%"
{res}{txt}file table1.xlsx saved

{com}. 
. putexcel H1 = "GCM"
{res}{txt}file table1.xlsx saved

{com}. putexcel H2 = "$\rho = 0.2$"
{res}{txt}file table1.xlsx saved

{com}. putexcel H3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel I3 = "%"
{res}{txt}file table1.xlsx saved

{com}. putexcel J2 = "$\rho = 0.5$"
{res}{txt}file table1.xlsx saved

{com}. putexcel J3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel K3 = "%"
{res}{txt}file table1.xlsx saved

{com}. putexcel L2 = "$\rho = 0.8$"
{res}{txt}file table1.xlsx saved

{com}. putexcel L3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel M3 = "%"
{res}{txt}file table1.xlsx saved

{com}. 
. putexcel N1 = "First Difference"
{res}{txt}file table1.xlsx saved

{com}. putexcel N2 = "$\rho = 0.2$"
{res}{txt}file table1.xlsx saved

{com}. putexcel N3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel O3 = "%"
{res}{txt}file table1.xlsx saved

{com}. putexcel P2 = "$\rho = 0.5$"
{res}{txt}file table1.xlsx saved

{com}. putexcel P3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel Q3 = "%"
{res}{txt}file table1.xlsx saved

{com}. putexcel R2 = "$\rho = 0.8$"
{res}{txt}file table1.xlsx saved

{com}. putexcel R3 = "coef"
{res}{txt}file table1.xlsx saved

{com}. putexcel S3 = "%"
{res}{txt}file table1.xlsx saved

{com}. 
. 
. ************************
. **T=5000; x1 = I(0), rho=.2, .5, .8
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *************************************************************************
. *************************************************************************
. *TABLE 1 - ADL
. *************************************************************************
. *************************************************************************
. set seed 457090451
{txt}
{com}. 
. **********
. **rho=.2
. **********
. *program drop combined_noq
. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg y l.y x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. 
. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}         97    .0246439    .0134346   .0001292   .0499567
{txt}
{com}. *_sim_1 = _b[L.y]
. *_sim_3 = _b[L.x1]
. sum _sim_1  _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000    .9989024    .0009467     .99029   1.000586
{txt}{space 6}_sim_3 {c |}{res}      2,000   -.9982733    .0141839  -1.055586  -.9482896
{txt}
{com}. 
. 
. *export results from ADL rho = 0.2 to excel file
. **$\hat{c -(}\alpha{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9989
{txt}
{com}. putexcel B4 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
{txt}
{com}. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}100
{txt}
{com}. putexcel C4 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9999
{txt}
{com}. putexcel B6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
{txt}
{com}. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel C6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_2$
. *ceof
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.9983
{txt}
{com}. putexcel B7 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
{txt}
{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel C7 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg y l.y x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}        111    .0257653    .0147279   .0002628    .049973
{txt}
{com}. *_sim_1 = _b[L.y]
. *_sim_3 = _b[L.x1]
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000    .9989108    .0009155   .9924613   1.000493
{txt}{space 7}_b_x1 {c |}{res}      2,000    .9999789    .0137387   .9546268   1.050255
{txt}{space 6}_sim_3 {c |}{res}      2,000   -.9994141    .0144915  -1.050412   -.952254
{txt}
{com}. 
. 
. *export results from ADL rho = 0.5 to excel file
. **$\hat{c -(}\alpha{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9989
{txt}
{com}. putexcel D4 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
{txt}
{com}. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}100
{txt}
{com}. putexcel E4 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0000
{txt}
{com}. putexcel D6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
{txt}
{com}. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel E6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_2$
. *coef
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.9994
{txt}
{com}. putexcel D7 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
{txt}
{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel E7 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg y l.y x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. 
. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}         97    .0236368    .0152229   .0001463   .0496133
{txt}
{com}. *_sim_1 = _b[L.y]
. *_sim_3 = _b[L.x1]
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000    .9988842    .0009084   .9939612   1.000751
{txt}{space 7}_b_x1 {c |}{res}      2,000    1.000413    .0142537   .9528446   1.057753
{txt}{space 6}_sim_3 {c |}{res}      2,000   -.9994568    .0144467  -1.044442  -.9584477
{txt}
{com}. 
. 
. 
. *export results from ADL rho = 0.8 to excel file
. **$\hat{c -(}\alpha{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9989
{txt}
{com}. putexcel F4 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_ly = abs(_sim_1/_sim_5) if abs(_sim_1/_sim_5)>=1.96
{txt}
{com}. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}100
{txt}
{com}. putexcel G4 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0004
{txt}
{com}. putexcel F6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1) >=1.96
{txt}
{com}. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel G6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_2$
. *coef
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.9995
{txt}
{com}. putexcel F7 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
{txt}
{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel G7 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. 
. *************************************************************************
. *************************************************************************
. *TABLE 1 - GECM
. *************************************************************************
. *************************************************************************
. *reset seed
. set seed 457090451
{txt}
{com}. 
. **********
. **rho=.2
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg d.y l.y d.x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}         97    .0246439    .0134346   .0001292   .0499567
{txt}
{com}. *_sim_1 = _b[L.y]
. *_sim_2 = _b[D.x1]
. *_sim_3 = _b[L.x1]
. sum  _sim_1 _sim_2 _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000   -.0010976    .0009467    -.00971   .0005858
{txt}{space 6}_sim_2 {c |}{res}      2,000    .9998638    .0140868   .9487547   1.048306
{txt}{space 6}_sim_3 {c |}{res}      2,000    .0015904    .0181709  -.0575766    .069093
{txt}
{com}. 
. 
. *export results from GECM rho = 0.2 to excel file
. **$\hat{c -(}\alpha^*{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.0011
{txt}
{com}. putexcel H5 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_ly = abs(_sim_1/_sim_6) if abs(_sim_1/_sim_6)>=1.96
{txt}(2,000 missing values generated)

{com}. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}0
{txt}
{com}. putexcel I5 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_1$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9999
{txt}
{com}. putexcel H6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_dx1 = abs(_sim_2/_sim_6) if abs(_sim_2/_sim_6) >=1.96
{txt}
{com}. quietly sum tstat_dx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel I6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta^*{c )-}_2$
. *coef
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}0.0016
{txt}
{com}. putexcel H8 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
{txt}(1,889 missing values generated)

{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}5.55
{txt}
{com}. putexcel I8 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg d.y l.y d.x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. 
. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}        111    .0257653    .0147279   .0002628    .049973
{txt}
{com}. *_sim_1 = _b[L.y]
. *_sim_2 = _b[D.x1]
. *_sim_3 = _b[L.x1]
. sum  _sim_1 _sim_2 _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000   -.0010892    .0009155  -.0075387   .0004932
{txt}{space 6}_sim_2 {c |}{res}      2,000    .9999789    .0137387   .9546268   1.050255
{txt}{space 6}_sim_3 {c |}{res}      2,000    .0005647    .0140469  -.0455149   .0422816
{txt}
{com}. 
. 
. 
. *export results from GECM rho = 0.5 to excel file
. **$\hat{c -(}\alpha^*{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.0011
{txt}
{com}. putexcel J5 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_ly = abs(_sim_1/_sim_6) if abs(_sim_1/_sim_6)>=1.96
{txt}(2,000 missing values generated)

{com}. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}0
{txt}
{com}. putexcel K5 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_1$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0000
{txt}
{com}. putexcel J6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_dx1 = abs(_sim_2/_sim_6) if abs(_sim_2/_sim_6) >=1.96
{txt}
{com}. quietly sum tstat_dx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel K6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta^*{c )-}_2$
. *coef
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}0.0006
{txt}
{com}. putexcel J8 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
{txt}(1,894 missing values generated)

{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}5.3
{txt}
{com}. putexcel K8 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg d.y l.y d.x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. 
. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}         97    .0236368    .0152229   .0001463   .0496133
{txt}
{com}. *_sim_1 = _b[L.y]
. *_sim_2 = _b[D.x1]
. *_sim_3 = _b[L.x1]
. sum  _sim_1 _sim_2 _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000   -.0011158    .0009084  -.0060388   .0007511
{txt}{space 6}_sim_2 {c |}{res}      2,000    1.000413    .0142537   .9528446   1.057753
{txt}{space 6}_sim_3 {c |}{res}      2,000    .0009564    .0089162  -.0281321   .0386523
{txt}
{com}. 
. 
. 
. *export results from GECM rho = 0.8 to excel file
. **$\hat{c -(}\alpha^*{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.0011
{txt}
{com}. putexcel L5 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_ly = abs(_sim_1/_sim_6) if abs(_sim_1/_sim_6)>=1.96
{txt}(2,000 missing values generated)

{com}. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}0
{txt}
{com}. putexcel M5 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta{c )-}_1$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0004
{txt}
{com}. putexcel L6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_dx1 = abs(_sim_2/_sim_6) if abs(_sim_2/_sim_6) >=1.96
{txt}
{com}. quietly sum tstat_dx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel M6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta^*{c )-}_2$
. *coef
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}0.0010
{txt}
{com}. putexcel L8 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_3/_sim_7) if abs(_sim_3/_sim_7)>=1.96
{txt}(1,897 missing values generated)

{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}5.15
{txt}
{com}. putexcel M8 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. *****************************************************************************
. *************************************************************************
. *TABLE 1 - FIRST-DIFFERENCE DV
. *************************************************************************
. *************************************************************************
. *reset seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.2
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg d.y x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}         96    .0253108    .0140307   .0001122   .0498891
{txt}
{com}. *_sim_2 = _b[L.x1]
. sum _b_x1 _sim_2

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      2,000    .9998685    .0140735   .9485242    1.04779
{txt}{space 6}_sim_2 {c |}{res}      2,000   -.9993665    .0141374  -1.056754  -.9488865
{txt}
{com}. 
. 
. *export results from First Difference rho = 0.2 to excel file
. **$\hat{c -(}\beta^{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9999
{txt}
{com}. putexcel N6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1)>=1.96
{txt}
{com}. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}100
{txt}
{com}. putexcel O6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta^{c )-}_2$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.9994
{txt}
{com}. putexcel N7 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_2/_sim_5) if abs(_sim_2/_sim_5) >=1.96
{txt}
{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel O7 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. 
. **********
. **rho=.5
. **********
. 
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg d.y x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}        107    .0252964    .0143134     .00027   .0492865
{txt}
{com}. *_sim_2 = _b[L.x1]
. sum _b_x1 _sim_2

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      2,000    .9999822    .0137413   .9547024   1.050099
{txt}{space 6}_sim_2 {c |}{res}      2,000     -1.0005    .0144768  -1.050957   -.953034
{txt}
{com}. 
. 
. *export results from First Difference rho = 0.5 to excel file
. **$\hat{c -(}\beta^{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0000
{txt}
{com}. putexcel P6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1)>=1.96
{txt}
{com}. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}100
{txt}
{com}. putexcel Q6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta^{c )-}_2$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}-1.0005
{txt}
{com}. putexcel P7 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_2/_sim_5) if abs(_sim_2/_sim_5) >=1.96
{txt}
{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel Q7 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. tsset t
{txt}  5{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  6{com}. gen x1=e1 if t==1
{txt}  7{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  8{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  9{com}. gen x2=u if t==1
{txt} 10{com}. replace x2=x2[_n-1] + u if t>1
{txt} 11{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 12{com}. gen y = x1 + x2
{txt} 13{com}. reg d.y x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(2000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_pvalu~g {c |}{res}        102    .0244508    .0154937   .0001391    .049013
{txt}
{com}. *_sim_2 = _b[L.x1]
. sum _b_x1 _sim_2

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      2,000    1.000419    .0142507   .9526367   1.057714
{txt}{space 6}_sim_2 {c |}{res}      2,000   -1.000566    .0144294  -1.045429   -.958738
{txt}
{com}. 
. 
. *export results from First Difference rho = 0.8 to excel file
. **$\hat{c -(}\beta^{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0004
{txt}
{com}. putexcel R6 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_x1 = abs(_b_x1/_se_x1) if abs(_b_x1/_se_x1)>=1.96
{txt}
{com}. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc' 
{res}100
{txt}
{com}. putexcel S6 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. **$\hat{c -(}\beta^{c )-}_2$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}-1.0006
{txt}
{com}. putexcel R7 = `r(mean)', nformat("#0.####")
{res}{txt}file table1.xlsx saved

{com}. *%
. gen tstat_lx1 = abs(_sim_2/_sim_5) if abs(_sim_2/_sim_5) >=1.96
{txt}
{com}. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel S7 = `perc'
{res}{txt}file table1.xlsx saved

{com}. 
. 
. 
. *************************************************************************
. *************************************************************************
. ****************************TABLE 2**************************************
. *************************************************************************
. *************************************************************************
. 
. **********************************
. *Table 2 layout
. **********************************
. *generate excel file to store Table 2 output
. putexcel set table2.xlsx, sheet(sheet1) replace
{res}{txt}Note: file will be replaced when the first {cmd:putexcel} command is issued

{com}. *add labels to table
. putexcel A1 = "$Y I(1), x^s_t I(0)$"
{res}{txt}file table2.xlsx saved

{com}. putexcel A2 = "$Y I(0), x^s_t I(0)$"
{res}{txt}file table2.xlsx saved

{com}. putexcel A3 = "$Y I(1), x^s_t I(1)$"
{res}{txt}file table2.xlsx saved

{com}. putexcel A4 = "$Y I(0), x^s_t I(1)$"
{res}{txt}file table2.xlsx saved

{com}. 
. 
. ************************
. **T=100; x1 = I(0), rho=.5
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *************************************************************************
. *************************************************************************
. *TABLE 2 - ADL
. *************************************************************************
. *************************************************************************
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 100
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.                 reg d.y x1 l.x1
{txt} 47{com}.                 mat P = r(p)
{txt} 48{com}.                 return scalar pvalue_bg = P[1,1] 
{txt} 49{com}.         
.         end
{txt}
{com}. 
.         
. simul _b _se p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 13 23 27 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 2 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 1 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}.                                 
. 
. 
. sum _b_x1 _sim_2

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      2,000    .9998595    .1013624   .6794183   1.365041
{txt}{space 6}_sim_2 {c |}{res}      2,000   -1.003063    .1030367  -1.334292  -.5862845
{txt}
{com}. gen tstat_x1_ = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1_ if tstat_x1_>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_x1_ {c |}{res}      2,000    9.908762    1.420192   5.333049   15.14183
{txt}
{com}. gen tstat_lx1 = abs(_sim_2/_sim_5)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_lx1 {c |}{res}      2,000    9.927531    1.447095   5.417834   15.21943
{txt}
{com}. sum _eq2_p_value_df_y if _eq2_p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}      1,571    .3897077    .2545861   .0501288   .9863328
{txt}
{com}. sum _eq2_p_value_df_y if _eq2_p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}        429    .0172043     .014755   .0000174   .0499682
{txt}
{com}. gen yI1=0
{txt}
{com}. replace yI1=1 if _eq2_p_value_df_y>.05
{txt}(1,571 real changes made)

{com}. local fsig_df_y = r(N)/10  
{txt}
{com}. sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}          7    .1476082    .1008611    .064952   .3593745
{txt}
{com}. sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}      1,993    .0002272    .0023005   2.95e-13   .0440806
{txt}
{com}. local fsig_df_x1 = r(N)/10 
{txt}
{com}. gen xI1=0
{txt}
{com}. replace xI1=1 if _eq2_p_value_df_x1>.05
{txt}(7 real changes made)

{com}.                                 
.                                 
. *Produce Table 2 crosstab of Y and X1 characteristics.
. tab yI1 xI1

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       428          1 {txt}{c |}{res}       429 
{txt}         1 {c |}{res}     1,565          6 {txt}{c |}{res}     1,571 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,993          7 {txt}{c |}{res}     2,000 

{txt}
{com}. gen y1x1=0
{txt}
{com}. replace y1x1 =1 if yI1==1 & xI1==1
{txt}(6 real changes made)

{com}. sum _b_x1 _sim_2 if y1x1==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}          6    1.054612    .0878624   .9428765   1.182716
{txt}{space 6}_sim_2 {c |}{res}          6   -1.020837    .0925836  -1.187428  -.9092814
{txt}
{com}. gen y1x0=0
{txt}
{com}. replace y1x0=1 if yI1==1 & xI1==0
{txt}(1,565 real changes made)

{com}. sum _b_x1 _sim_2 if y1x0==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      1,565    .9961916    .1016562   .6794183   1.365041
{txt}{space 6}_sim_2 {c |}{res}      1,565   -.9997917    .1021844  -1.332114  -.5862845
{txt}
{com}. gen tstat_x1__ = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1__ if tstat_x1_>=1.96               

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 2}tstat_x1__ {c |}{res}      2,000    9.908762    1.420192   5.333049   15.14183
{txt}
{com}.                                 
. *Turn crosstabs into matrix to export Y and X1 characteristics into excel                               
. tab yI1 xI1, matcell(x)                 

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       428          1 {txt}{c |}{res}       429 
{txt}         1 {c |}{res}     1,565          6 {txt}{c |}{res}     1,571 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,993          7 {txt}{c |}{res}     2,000 

{txt}
{com}. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  428     1
{txt}r2 {res} 1565     6
{reset}
{com}. 
. scalar a = x[2,1]
{txt}
{com}. local a = a/2000
{txt}
{com}. di `a'
{res}.7825
{txt}
{com}. putexcel B1 = `a'
{res}{txt}file table2.xlsx saved

{com}. 
. scalar b = x[1,1]
{txt}
{com}. local b = b/2000
{txt}
{com}. di `b'
{res}.214
{txt}
{com}. putexcel B2 = `b'
{res}{txt}file table2.xlsx saved

{com}. 
. scalar c = x[2,2]
{txt}
{com}. local c = c/2000
{txt}
{com}. di `c'
{res}.003
{txt}
{com}. putexcel B3 = `c'
{res}{txt}file table2.xlsx saved

{com}. 
. scalar d = x[1,2]
{txt}
{com}. local d = d/2000
{txt}
{com}. di `d'
{res}.0005
{txt}
{com}. putexcel B4 = `d'
{res}{txt}file table2.xlsx saved

{com}. 
. 
. 
. *************************************************************************
. *************************************************************************
. ****************************TABLE 3**************************************
. *************************************************************************
. *************************************************************************
. 
. 
. **********************************
. *Table 3 layout
. **********************************
. *generate excel file to store Table 3 output
. putexcel set table3.xlsx, sheet(sheet1) replace
{res}{txt}Note: file will be replaced when the first {cmd:putexcel} command is issued

{com}. *add labels to table
. 
. putexcel B1 = "First Difference, Y = I(1), $x^s_t$ = I (0)"
{res}{txt}file table3.xlsx saved

{com}. putexcel B2 = "coef"
{res}{txt}file table3.xlsx saved

{com}. putexcel C2 = "%"
{res}{txt}file table3.xlsx saved

{com}. putexcel D1 = "ADL, Y = I(0), $x^s_t$ = I(0)"
{res}{txt}file table3.xlsx saved

{com}. putexcel D2 = "coef"
{res}{txt}file table3.xlsx saved

{com}. putexcel E2 = "%"
{res}{txt}file table3.xlsx saved

{com}. putexcel A3 = "$\hat{c -(}\alpha{c )-}_1$"
{res}{txt}file table3.xlsx saved

{com}. putexcel A4 = "$\hat{c -(}\beta{c )-}_1$"
{res}{txt}file table3.xlsx saved

{com}. putexcel A5 = "$\hat{c -(}\beta{c )-}_2$"
{res}{txt}file table3.xlsx saved

{com}. putexcel A6 = "% cases"
{res}{txt}file table3.xlsx saved

{com}. 
. 
. *Estimation of relationship between DV and IV where former is I(1) and latter is I(0)* 
. 
. ************************
. **T=100; x1 = I(0), rho=.5
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *************************************************************************
. *************************************************************************
. *Differenced DV
. *************************************************************************
. *************************************************************************
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 100
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.                 reg d.y x1 l.x1
{txt} 47{com}.                 mat P = r(p)
{txt} 48{com}.                 return scalar pvalue_bg = P[1,1] 
{txt} 49{com}.                 
.         end
{txt}
{com}. 
.         
.                                 simul _b _se p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 13 23 27 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 2 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 1 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 sum _b_x1 _sim_2

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      2,000    .9998595    .1013624   .6794183   1.365041
{txt}{space 6}_sim_2 {c |}{res}      2,000   -1.003063    .1030367  -1.334292  -.5862845
{txt}
{com}.                                 gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}.                                 sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 4}tstat_x1 {c |}{res}      2,000    9.908762    1.420192   5.333049   15.14183
{txt}
{com}.                                 gen tstat_lx1 = abs(_sim_2/_sim_5)
{txt}
{com}.                                 sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_lx1 {c |}{res}      2,000    9.927531    1.447095   5.417834   15.21943
{txt}
{com}.                                 sum _eq2_p_value_df_y if _eq2_p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}      1,571    .3897077    .2545861   .0501288   .9863328
{txt}
{com}.                                 sum _eq2_p_value_df_y if _eq2_p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}        429    .0172043     .014755   .0000174   .0499682
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if _eq2_p_value_df_y>.05
{txt}(1,571 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}          7    .1476082    .1008611    .064952   .3593745
{txt}
{com}.                                 sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}      1,993    .0002272    .0023005   2.95e-13   .0440806
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if _eq2_p_value_df_x1>.05
{txt}(7 real changes made)

{com}.                                 tab yI1 xI1

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       428          1 {txt}{c |}{res}       429 
{txt}         1 {c |}{res}     1,565          6 {txt}{c |}{res}     1,571 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,993          7 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 gen y1x1=0
{txt}
{com}.                                 replace y1x1 =1 if yI1==1 & xI1==1
{txt}(6 real changes made)

{com}.                                 sum _b_x1 _sim_2 if y1x1==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}          6    1.054612    .0878624   .9428765   1.182716
{txt}{space 6}_sim_2 {c |}{res}          6   -1.020837    .0925836  -1.187428  -.9092814
{txt}
{com}.                                 gen y1x0=0
{txt}
{com}.                                 replace y1x0=1 if yI1==1 & xI1==0
{txt}(1,565 real changes made)

{com}.                                 sum _b_x1 _sim_2 if y1x0==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      1,565    .9961916    .1016562   .6794183   1.365041
{txt}{space 6}_sim_2 {c |}{res}      1,565   -.9997917    .1021844  -1.332114  -.5862845
{txt}
{com}. 
.                                 sum tstat_x1 if y1x0==1                 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 4}tstat_x1 {c |}{res}      1,565    9.846167    1.410536   5.333049   15.14183
{txt}
{com}.                                 sum tstat_lx1 if y1x0==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_lx1 {c |}{res}      1,565    9.868092    1.421935   5.417834   15.21943
{txt}
{com}. 
. *Mean estimate of IV coefficients across simulations and frequency of statistical significance*
.                                 sum _b_x1 _sim_2 if y1x0==1 & tstat_x1>1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}      1,565    .9961916    .1016562   .6794183   1.365041
{txt}{space 6}_sim_2 {c |}{res}      1,565   -.9997917    .1021844  -1.332114  -.5862845
{txt}
{com}. 
. 
. 
.                                 
. *export results from First Differenced DV rho = 0.5 to excel file
. 
. **$\hat{c -(}\beta^{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9999
{txt}
{com}. putexcel B4 = `r(mean)', nformat("#0.####")
{res}{txt}file table3.xlsx saved

{com}. *%
. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel C4 = `perc'
{res}{txt}file table3.xlsx saved

{com}. 
. 
.                                 
. **$\hat{c -(}\beta^{c )-}_2$
. *coef
. quietly sum _sim_2
{txt}
{com}. di %6.4f `r(mean)'
{res}-1.0031
{txt}
{com}. putexcel B5 = `r(mean)', nformat("#0.####")
{res}{txt}file table3.xlsx saved

{com}. *%
. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel C5 = `perc'
{res}{txt}file table3.xlsx saved

{com}. 
. 
. 
. *Estimation of relationship between DV and IV where former is I(0) and latter is I(0)* 
. 
. ************************
. **T=100; x1 = I(0), rho=.5
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. 
. *************************************************************************
. *************************************************************************
. *ADL
. *************************************************************************
. *************************************************************************
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 100
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.                 reg y l.y x1 l.x1
{txt} 47{com}.                 mat P = r(p)
{txt} 48{com}.                 return scalar pvalue_bg = P[1,1] 
{txt} 49{com}.                 
.         end
{txt}
{com}. 
.         
.                                 simul _b _se p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 13 23 27 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 2 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 1 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}      2,000    .9460526    .0446224   .6919978   1.029792
{txt}{space 7}_b_x1 {c |}{res}      2,000     1.00008    .1011249   .6756603   1.367891
{txt}{space 6}_sim_3 {c |}{res}      2,000   -.9490748    .1126602  -1.299508  -.5486233
{txt}
{com}.                                 gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}.                                 sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 4}tstat_ly {c |}{res}      2,000    35.58593    15.14379   9.404983   108.8567
{txt}
{com}.                                 gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}.                                 sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 4}tstat_x1 {c |}{res}      2,000    9.963066    1.433587   5.502026   15.64954
{txt}
{com}.                                 gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}.                                 sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_lx1 {c |}{res}      2,000     8.99611    1.485481   4.018018   15.25006
{txt}
{com}.                                 sum _eq2_p_value_df_y if _eq2_p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}      1,571    .3897077    .2545861   .0501288   .9863328
{txt}
{com}.                                 sum _eq2_p_value_df_y if _eq2_p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}        429    .0172043     .014755   .0000174   .0499682
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if _eq2_p_value_df_y>.05
{txt}(1,571 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}          7    .1476082    .1008611    .064952   .3593745
{txt}
{com}.                                 sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}      1,993    .0002272    .0023005   2.95e-13   .0440806
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if _eq2_p_value_df_x1>.05
{txt}(7 real changes made)

{com}.                                 tab yI1 xI1

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       428          1 {txt}{c |}{res}       429 
{txt}         1 {c |}{res}     1,565          6 {txt}{c |}{res}     1,571 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,993          7 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 gen y1x1=0
{txt}
{com}.                                 replace y1x1 =1 if yI1==1 & xI1==1
{txt}(6 real changes made)

{com}.                                 sum _b_x1 _sim_3 if y1x1==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 7}_b_x1 {c |}{res}          6    1.045441    .0903427    .947088   1.167078
{txt}{space 6}_sim_3 {c |}{res}          6   -.9941909     .103262  -1.193589  -.8910931
{txt}
{com}.                                 gen y1x0=0
{txt}
{com}.                                 replace y1x0=1 if yI1==1 & xI1==0
{txt}(1,565 real changes made)

{com}.                                 gen y0x0=0
{txt}
{com}.                                 replace y0x0=1 if yI1==0 & xI1==0
{txt}(428 real changes made)

{com}.                                 sum _sim_1 _b_x1 _sim_3 if y0x0==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}        428    .8836592    .0417294   .6919978   .9566926
{txt}{space 7}_b_x1 {c |}{res}        428    1.008312    .0981621    .691892   1.319694
{txt}{space 6}_sim_3 {c |}{res}        428   -.9011232    .1180325  -1.198269  -.5549219
{txt}
{com}.                                 sum tstat_ly if y0x0==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 4}tstat_ly {c |}{res}        428     20.7937    5.820769   9.404983   49.36953
{txt}
{com}.                                 sum tstat_x1 if y0x0==1                 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 4}tstat_x1 {c |}{res}        428    10.32232    1.461948   6.207458   14.76059
{txt}
{com}.                                 sum tstat_lx1 if y0x0==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_lx1 {c |}{res}        428    8.431383    1.540679   4.018018   14.02993
{txt}
{com}.                                 
. *Mean estimate of IV coefficients across simulations and frequency of statistical significance*
.                                 sum _sim_1 _b_x1 _sim_3 if y0x0==1 & tstat_x1>1.96      

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_sim_1 {c |}{res}        428    .8836592    .0417294   .6919978   .9566926
{txt}{space 7}_b_x1 {c |}{res}        428    1.008312    .0981621    .691892   1.319694
{txt}{space 6}_sim_3 {c |}{res}        428   -.9011232    .1180325  -1.198269  -.5549219
{txt}
{com}.                                 
. 
.                                 
. *export results from ADL rho = 0.5 to excel file
. 
. **$\hat{c -(}\alpha^{c )-}_1$
. *coef
. quietly sum _sim_1
{txt}
{com}. di %6.4f `r(mean)'
{res}0.9461
{txt}
{com}. putexcel D3 = `r(mean)', nformat("#0.####")
{res}{txt}file table3.xlsx saved

{com}. *%
. quietly sum tstat_ly
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel E3 = `perc'
{res}{txt}file table3.xlsx saved

{com}.                                 
. **$\hat{c -(}\beta^{c )-}_1$
. *coef
. quietly sum _b_x1
{txt}
{com}. di %6.4f `r(mean)'
{res}1.0001
{txt}
{com}. putexcel D4 = `r(mean)', nformat("#0.####")
{res}{txt}file table3.xlsx saved

{com}. *%
. quietly sum tstat_x1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel E4 = `perc'
{res}{txt}file table3.xlsx saved

{com}. 
.                                 
. **$\hat{c -(}\beta^{c )-}_2$
. *coef
. quietly sum _sim_3
{txt}
{com}. di %6.4f `r(mean)'
{res}-0.9491
{txt}
{com}. putexcel D5 = `r(mean)', nformat("#0.####")
{res}{txt}file table3.xlsx saved

{com}. *%
. quietly sum tstat_lx1
{txt}
{com}. local perc = 100*`r(N)'/2000
{txt}
{com}. di `perc'
{res}100
{txt}
{com}. putexcel E5 = `perc'
{res}{txt}file table3.xlsx saved

{com}. 
. 
. *export the percentages that correctly identify the series to excel
. 
. *Turn crosstabs into matrix to export Y and X1 characteristics into excel                               
. tab yI1 xI1, matcell(x)                 

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       428          1 {txt}{c |}{res}       429 
{txt}         1 {c |}{res}     1,565          6 {txt}{c |}{res}     1,571 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,993          7 {txt}{c |}{res}     2,000 

{txt}
{com}. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  428     1
{txt}r2 {res} 1565     6
{reset}
{com}. 
. scalar a = x[2,1]
{txt}
{com}. local a = a/2000
{txt}
{com}. di `a'
{res}.7825
{txt}
{com}. putexcel C6 = `a'
{res}{txt}file table3.xlsx saved

{com}. 
. scalar b = x[1,1]
{txt}
{com}. local b = b/2000
{txt}
{com}. di `b'
{res}.214
{txt}
{com}. putexcel E6 = `b'
{res}{txt}file table3.xlsx saved

{com}. 
. 
. *************************************************************************
. *************************************************************************
. ****************************TABLE 4**************************************
. *************************************************************************
. *************************************************************************
. 
. ***************
. *Table 4 layout
. ***************
. 
. **********************************
. *generate excel file to store Table 4 output
. putexcel set table4.xlsx, sheet(sheet1) replace
{res}{txt}Note: file will be replaced when the first {cmd:putexcel} command is issued

{com}. *add labels to table
. *row labels
. putexcel A3 = "Y I(1), X I(0)"
{res}{txt}file table4.xlsx saved

{com}. putexcel A4 = "Y I(0), X I(0)"
{res}{txt}file table4.xlsx saved

{com}. putexcel A5 = "Y I(1), X I(1)"
{res}{txt}file table4.xlsx saved

{com}. putexcel A6 = "Y I(0), X I(1)"
{res}{txt}file table4.xlsx saved

{com}. 
. *column labels
. putexcel B1 = "T = 50"
{res}{txt}file table4.xlsx saved

{com}. putexcel B2 = "$\rho = 0.2$"
{res}{txt}file table4.xlsx saved

{com}. putexcel C2 = "$\rho = 0.5$"
{res}{txt}file table4.xlsx saved

{com}. putexcel D2 = "$\rho = 0.8$"
{res}{txt}file table4.xlsx saved

{com}. 
. putexcel E1 = "T = 100"
{res}{txt}file table4.xlsx saved

{com}. putexcel E2 = "$\rho = 0.2$"
{res}{txt}file table4.xlsx saved

{com}. putexcel F2 = "$\rho = 0.5$"
{res}{txt}file table4.xlsx saved

{com}. putexcel G2 = "$\rho = 0.8$"
{res}{txt}file table4.xlsx saved

{com}. 
. putexcel H1 = "T = 200"
{res}{txt}file table4.xlsx saved

{com}. putexcel H2 = "$\rho = 0.2$"
{res}{txt}file table4.xlsx saved

{com}. putexcel I2 = "$\rho = 0.5$"
{res}{txt}file table4.xlsx saved

{com}. putexcel J2 = "$\rho = 0.8$"
{res}{txt}file table4.xlsx saved

{com}. 
. ************************
. **T=50; x1 = I(0), rho=.2
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. **********
. **rho=.2
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 50
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,520    .3868172     .259913   .0502341    .998079
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        480    .0169995    .0151128   4.81e-07   .0497743
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,520 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}          9    .1169302    .0635736   .0595643    .245712
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,991    .0003179    .0025506   1.18e-16   .0481139
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(9 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       480          0 {txt}{c |}{res}       480 
{txt}         1 {c |}{res}     1,511          9 {txt}{c |}{res}     1,520 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,991          9 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
. 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  480     0
{txt}r2 {res} 1511     9
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}75.55
{txt}
{com}. putexcel B3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di b
{res}24
{txt}
{com}. putexcel B4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}.45
{txt}
{com}. putexcel B5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}0
{txt}
{com}. putexcel B6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. 
. 
. ************************
. **T=100; x1 = I(0), rho=.2
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. **********
. **rho=.2
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 100
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,475    .3859048    .2605866   .0505452   .9821792
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        525    .0149723    .0151973   1.93e-07   .0498044
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,475 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}          3    .1159335    .0144774   .0992209   .1246253
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,997    .0001257    .0018508   4.12e-21   .0451945
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(3 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                         
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       525          0 {txt}{c |}{res}       525 
{txt}         1 {c |}{res}     1,472          3 {txt}{c |}{res}     1,475 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,997          3 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
. 
. 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  525     0
{txt}r2 {res} 1472     3
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}73.6
{txt}
{com}. putexcel E3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}73.6
{txt}
{com}. putexcel E4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}.15
{txt}
{com}. putexcel E5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}0
{txt}
{com}. putexcel E6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. ************************
. **T=200; x1 = I(0), rho=.2
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.2
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 200
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,528    .3980878    .2629234   .0500251   .9971423
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        472    .0151127    .0150904   9.55e-07   .0499554
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,528 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}          1    .0619604           .   .0619604   .0619604
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,999    .0000185    .0006518   1.29e-27   .0287234
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(1 real change made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell (x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       472          0 {txt}{c |}{res}       472 
{txt}         1 {c |}{res}     1,527          1 {txt}{c |}{res}     1,528 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,999          1 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
.                                 
.                                 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  472     0
{txt}r2 {res} 1527     1
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}76.35
{txt}
{com}. putexcel H3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}76.35
{txt}
{com}. putexcel H4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}.05
{txt}
{com}. putexcel H5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}0
{txt}
{com}. putexcel H6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. ************************
. **T=50; x1 = I(0), rho=.5
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 50
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,660    .4018657    .2568455   .0500421   .9964736
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        340    .0201478    .0152614   4.75e-09   .0497231
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,660 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}         73    .1183534    .0860446   .0512919   .3886483
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,927    .0037584    .0070392   2.08e-10   .0495206
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(73 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       330         10 {txt}{c |}{res}       340 
{txt}         1 {c |}{res}     1,597         63 {txt}{c |}{res}     1,660 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,927         73 {txt}{c |}{res}     2,000 

{txt}
{com}. 
.                                 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
.                                 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  330    10
{txt}r2 {res} 1597    63
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}79.85
{txt}
{com}. putexcel C3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}79.85
{txt}
{com}. putexcel C4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}3.15
{txt}
{com}. putexcel C5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}.5
{txt}
{com}. putexcel C6 = d
{res}{txt}file table4.xlsx saved

{com}.                         
. 
. ************************
. **T=100; x1 = I(0), rho=.5
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 100
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,571    .3897077    .2545861   .0501288   .9863328
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        429    .0172043     .014755   .0000174   .0499682
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,571 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}          7    .1476082    .1008611    .064952   .3593745
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,993    .0002272    .0023005   2.95e-13   .0440806
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(7 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       428          1 {txt}{c |}{res}       429 
{txt}         1 {c |}{res}     1,565          6 {txt}{c |}{res}     1,571 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,993          7 {txt}{c |}{res}     2,000 

{txt}
{com}. 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
.                                 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  428     1
{txt}r2 {res} 1565     6
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}78.25
{txt}
{com}. putexcel F3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}78.25
{txt}
{com}. putexcel F4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}.3
{txt}
{com}. putexcel F5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}.05
{txt}
{com}. putexcel F6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. ************************
. **T=200; x1 = I(0), rho=.5
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.5
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 200
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,549    .3840419    .2567653   .0503551   .9902613
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        451    .0158758    .0146532   4.73e-06   .0497479
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,549 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}          1     .094486           .    .094486    .094486
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,999    .0000206    .0007117   2.79e-20    .031644
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(1 real change made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell (x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       451          0 {txt}{c |}{res}       451 
{txt}         1 {c |}{res}     1,548          1 {txt}{c |}{res}     1,549 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,999          1 {txt}{c |}{res}     2,000 

{txt}
{com}. 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
. 
.                                 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  451     0
{txt}r2 {res} 1548     1
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}77.4
{txt}
{com}. putexcel I3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}77.4
{txt}
{com}. putexcel I4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}.05
{txt}
{com}. putexcel I5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}0
{txt}
{com}. putexcel I6 = d
{res}{txt}file table4.xlsx saved

{com}. 
.         
. 
. ************************
. **T=50; x1 = I(0), rho=.8
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 50
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,827    .4424665    .2600039   .0503814   .9968744
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        173    .0219914    .0144353   .0000549   .0499028
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,827 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,332    .2069281    .1522136   .0500149   .9372349
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}        668    .0192206    .0149651   .0000144   .0499664
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(1,332 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}        74         99 {txt}{c |}{res}       173 
{txt}         1 {c |}{res}       594      1,233 {txt}{c |}{res}     1,827 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}       668      1,332 {txt}{c |}{res}     2,000 

{txt}
{com}. 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
.                                 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}   74    99
{txt}r2 {res}  594  1233
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}29.7
{txt}
{com}. putexcel D3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di b
{res}3.7
{txt}
{com}. putexcel D4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}61.65
{txt}
{com}. putexcel D5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}4.95
{txt}
{com}. putexcel D6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. *using critical value of 0.1 (footnote 14)
. 
. 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.1 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,668    .4776491     .244537   .1012668   .9968744
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        332    .0466027    .0294873   .0000549   .0996849
{txt}
{com}.                                 gen yI1_=0
{txt}
{com}.                                 replace yI1_=1 if p_value_df_y>.1
{txt}(1,668 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}        946    .2616426    .1490287   .1001754   .9372349
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,054    .0388555    .0297159   .0000144   .0997807
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1_=0
{txt}
{com}.                                 replace xI1_=1 if p_value_df_x1>.1
{txt}(946 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1_ xI1_, matcell(x)

           {txt}{c |}         xI1_
      yI1_ {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       197        135 {txt}{c |}{res}       332 
{txt}         1 {c |}{res}       857        811 {txt}{c |}{res}     1,668 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,054        946 {txt}{c |}{res}     2,000 

{txt}
{com}. 
.                                 
. *Turn crosstab into matrix 
.                                 
. matrix list x
{res}
{txt}x[2,2]
     c1   c2
r1 {res} 197  135
{txt}r2 {res} 857  811
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}42.85
{txt}
{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di b
{res}9.85
{txt}
{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}40.55
{txt}
{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}6.75
{txt}
{com}. 
. 
. 
. 
. 
. ************************
. **T=100; x1 = I(0), rho=.8
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 100
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,732    .4121211    .2561904   .0501311   .9912249
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        268    .0205925    .0143349   .0001054   .0498037
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,732 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}        296    .1099673    .0745814    .050225   .6435528
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,704    .0104721    .0119848   9.46e-07   .0498235
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(296 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       233         35 {txt}{c |}{res}       268 
{txt}         1 {c |}{res}     1,471        261 {txt}{c |}{res}     1,732 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,704        296 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
. 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  233    35
{txt}r2 {res} 1471   261
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}73.55
{txt}
{com}. putexcel G3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}73.55
{txt}
{com}. putexcel G4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}13.05
{txt}
{com}. putexcel G5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}1.75
{txt}
{com}. putexcel G6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. ************************
. **T=200; x1 = I(0), rho=.8
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *reset the seed
. set seed 457090451
{txt}
{com}. 
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 200
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  9{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 10{com}. gen x2=u if t==1
{txt} 11{com}. replace x2=x2[_n-1] + u if t>1
{txt} 12{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 13{com}. gen y = x1 + x2
{txt} 14{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 15{com}.                 predict resid_y, r 
{txt} 16{com}.                 wntestq resid_y 
{txt} 17{com}.                 local i=r(p)  
{txt} 18{com}.                 drop resid_y 
{txt} 19{com}.                 while `i'<0.05 {c -(}  
{txt} 20{com}.                 local j =`j'+ 1  
{txt} 21{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 22{com}.                 predict resid_y, r
{txt} 23{com}.                 wntestq resid_y 
{txt} 24{com}.                 local i=r(p)  
{txt} 25{com}.                 drop resid_y 
{txt} 26{com}.                 {c )-} 
{txt} 27{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 28{com}.                 predict resid_x1, r 
{txt} 29{com}.                 wntestq resid_x1 
{txt} 30{com}.                 local k=r(p)  
{txt} 31{com}.                 drop resid_x1 
{txt} 32{com}.                 while `k'<0.05 {c -(}  
{txt} 33{com}.                 local l=`l'+ 1  
{txt} 34{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 35{com}.                 predict resid_x1, r 
{txt} 36{com}.                 wntestq resid_x1 
{txt} 37{com}.                 local k=r(p)  
{txt} 38{com}.                 drop resid_x1 
{txt} 39{com}.                 {c )-} 
{txt} 40{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 41{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 42{com}.                 return scalar df_y=r(Zt)
{txt} 43{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 44{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 45{com}.                 return scalar df_x1=r(Zt)
{txt} 46{com}.                 
.         end
{txt}
{com}. 
.         
. 
.                                 simul p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 10 20 24 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 5 20 24 2}{...}
{p2col :p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 4 20 24 2}{...}
{p2col :p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 
. *Diagnose I(1) of DV*
.                                 sum p_value_df_y if p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}      1,659    .3946394    .2521406   .0502206   .9915859
{txt}
{com}.                                 sum p_value_df_y if p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df_y {c |}{res}        341    .0193122    .0152231   .0001771   .0498958
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if p_value_df_y>.05
{txt}(1,659 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 
. *Diagnose I(1) of IV*
.                                 sum p_value_df_x1 if p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}          4    .1589463    .1723121   .0592756   .4166673
{txt}
{com}.                                 sum p_value_df_x1 if p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
p_value_df~1 {c |}{res}      1,996    .0005167    .0021419   2.90e-10   .0408673
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if p_value_df_x1>.05
{txt}(4 real changes made)

{com}. 
. *Tabulate characteristics of DV and IV*
.                                 tab yI1 xI1, matcell(x)

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}       341          0 {txt}{c |}{res}       341 
{txt}         1 {c |}{res}     1,655          4 {txt}{c |}{res}     1,659 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}     1,996          4 {txt}{c |}{res}     2,000 

{txt}
{com}. 
. *Turn crosstab into matrix to export DV and IV characteristics into excel
. 
. matrix list x
{res}
{txt}x[2,2]
      c1    c2
r1 {res}  341     0
{txt}r2 {res} 1655     4
{reset}
{com}. 
. scalar a = x[2,1]/2000*100
{txt}
{com}. di a
{res}82.75
{txt}
{com}. putexcel J3 = a
{res}{txt}file table4.xlsx saved

{com}. 
. scalar b = x[1,1]/2000*100
{txt}
{com}. di a
{res}82.75
{txt}
{com}. putexcel J4 = b
{res}{txt}file table4.xlsx saved

{com}. 
. scalar c = x[2,2]/2000*100
{txt}
{com}. di c
{res}.2
{txt}
{com}. putexcel J5 = c
{res}{txt}file table4.xlsx saved

{com}. 
. scalar d = x[1,2]/2000*100
{txt}
{com}. di d
{res}0
{txt}
{com}. putexcel J6 = d
{res}{txt}file table4.xlsx saved

{com}. 
. 
. 
. 
. **************************************************************************************
. **************************************************************************************
. ****************************TABLE 4 - SUPPLEMENT**************************************
. **************************************************************************************
. **************************************************************************************
. 
. 
. 
. *Estimation of relationship between DV and IV where former seems I(1) and latter seems I(1)* 
. 
. ************************
. **T=50; x1 = I(0), rho=.8
. **x2 = I(1)
. **y = x1 + x2 
. *************************
. 
. *************************************************************************
. *************************************************************************
. *Differenced specification
. *************************************************************************
. *************************************************************************
. *keeping the same seed as used above
. 
. **********
. **rho=.8
. **********
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) r(real 1) s(real 0) v(real 0) w(real 0) obs(integer 5000)] 
{txt}  2{com}. *#d
. drop _all
{txt}  3{com}. set obs 50
{txt}  4{com}. gen t = _n
{txt}  5{com}. tsset t
{txt}  6{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  7{com}. gen x1=e1 if t==1
{txt}  8{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  9{com}. gen dx1=x1-l.x1
{txt} 10{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt} 11{com}. gen x2=u if t==1
{txt} 12{com}. replace x2=x2[_n-1] + u if t>1
{txt} 13{com}. *gen combined times series (z) that is a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 14{com}. gen y = x1 + x2
{txt} 15{com}. 
. 
.                 
.                 reg d.y l.y  
{txt} 16{com}.                 predict resid_y, r 
{txt} 17{com}.                 wntestq resid_y 
{txt} 18{com}.                 local i=r(p)  
{txt} 19{com}.                 drop resid_y 
{txt} 20{com}.                 while `i'<0.05 {c -(}  
{txt} 21{com}.                 local j =`j'+ 1  
{txt} 22{com}.                 reg d.y l.y l(1/`j')d.y  
{txt} 23{com}.                 predict resid_y, r
{txt} 24{com}.                 wntestq resid_y 
{txt} 25{com}.                 local i=r(p)  
{txt} 26{com}.                 drop resid_y 
{txt} 27{com}.                 {c )-} 
{txt} 28{com}.                 
.                 
.                 reg d.x1 l.x1  
{txt} 29{com}.                 predict resid_x1, r 
{txt} 30{com}.                 wntestq resid_x1 
{txt} 31{com}.                 local k=r(p)  
{txt} 32{com}.                 drop resid_x1 
{txt} 33{com}.                 while `k'<0.05 {c -(}  
{txt} 34{com}.                 local l=`l'+ 1  
{txt} 35{com}.                 reg d.x1 l.x1 l(1/`l')d.x1  
{txt} 36{com}.                 predict resid_x1, r 
{txt} 37{com}.                 wntestq resid_x1 
{txt} 38{com}.                 local k=r(p)  
{txt} 39{com}.                 drop resid_x1 
{txt} 40{com}.                 {c )-} 
{txt} 41{com}.                 
.         
.                 
.                 dfuller y, lags(`j')  
{txt} 42{com}.                 
.                 return scalar p_value_df_y=r(p)
{txt} 43{com}.                 return scalar df_y=r(Zt)
{txt} 44{com}.                 
.                 dfuller x1, lags(`l') 
{txt} 45{com}.                 return scalar p_value_df_x1=r(p) 
{txt} 46{com}.                 return scalar df_x1=r(Zt)
{txt} 47{com}.                 
.                 reg d.y dx1
{txt} 48{com}.                 mat P = r(p)
{txt} 49{com}.                 return scalar pvalue_bg = P[1,1] 
{txt} 50{com}.                 
.         end
{txt}
{com}. 
.         
.                                 simul _b _se p_value_df_y=r(p_value_df_y) p_value_df_x1=r(p_value_df_x1), reps(2000) nodots: combined_noq 
{p2colset 13 23 27 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 2 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_y:}{res:r(p_value_df_y)}{p_end}
{p2colset 1 23 27 2}{...}
{p2col :{txt}[{res:_eq2}]p_value_df_x1:}{res:r(p_value_df_x1)}{p_end}


{com}. 
.                                 sum _b_dx1 _se_dx1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_b_dx1 {c |}{res}      2,000    .9985427     .138145    .582692   1.501121
{txt}{space 5}_se_dx1 {c |}{res}      2,000    .1374755    .0202222   .0765309   .2360806
{txt}
{com}. 
.                                 sum _eq2_p_value_df_y if _eq2_p_value_df_y>0.05 

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}      1,808    .4402527    .2647384   .0500204   .9989168
{txt}
{com}.                                 sum _eq2_p_value_df_y if _eq2_p_value_df_y<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~y {c |}{res}        192    .0209585    .0147485   .0000354   .0496074
{txt}
{com}.                                 gen yI1=0
{txt}
{com}.                                 replace yI1=1 if _eq2_p_value_df_y>.05
{txt}(1,808 real changes made)

{com}.                                 local fsig_df_y = r(N)/10  
{txt}
{com}.                                 sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1>0.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}      1,266    .2080418    .1493648    .050155   .9577882
{txt}
{com}.                                 sum _eq2_p_value_df_x1 if _eq2_p_value_df_x1<.05

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
_eq2_p_val~1 {c |}{res}        734    .0191843    .0148936   .0000238   .0499204
{txt}
{com}.                                 local fsig_df_x1 = r(N)/10 
{txt}
{com}.                                 gen xI1=0
{txt}
{com}.                                 replace xI1=1 if _eq2_p_value_df_x1>.05
{txt}(1,266 real changes made)

{com}.                                 tab yI1 xI1

           {txt}{c |}          xI1
       yI1 {c |}         0          1 {c |}     Total
{hline 11}{c +}{hline 22}{c +}{hline 10}
         0 {c |}{res}        88        104 {txt}{c |}{res}       192 
{txt}         1 {c |}{res}       646      1,162 {txt}{c |}{res}     1,808 
{txt}{hline 11}{c +}{hline 22}{c +}{hline 10}
     Total {c |}{res}       734      1,266 {txt}{c |}{res}     2,000 

{txt}
{com}.                                 gen y1x0=0
{txt}
{com}.                                 replace y1x0=1 if yI1==1 & xI1==0
{txt}(646 real changes made)

{com}.                                 gen y1x1=0
{txt}
{com}.                                 replace y1x1 =1 if yI1==1 & xI1==1
{txt}(1,162 real changes made)

{com}.                                 sum _b_dx1 if y1x1==1

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_b_dx1 {c |}{res}      1,162    .9991382    .1393749    .582692   1.474369
{txt}
{com}.                                 gen tstat_dx1 = abs(_b_dx1/_se_dx1)
{txt}
{com}.                                 sum tstat_dx1 if tstat_dx1>=1.96        

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_dx1 {c |}{res}      2,000    7.425754    1.510623   2.589466   13.35536
{txt}
{com}.                                 sum tstat_dx1 if y1x1==1                

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 3}tstat_dx1 {c |}{res}      1,162    7.366139    1.505689   2.589466   13.18868
{txt}
{com}.                                 sum _b_dx1 tstat_dx1 if y1x1==1 & tstat_dx1>1.96

{txt}    Variable {c |}        Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 57}
{space 6}_b_dx1 {c |}{res}      1,162    .9991382    .1393749    .582692   1.474369
{txt}{space 3}tstat_dx1 {c |}{res}      1,162    7.366139    1.505689   2.589466   13.18868
{txt}
{com}.         
. 
. 
. 
{txt}end of do-file

{com}. log close
      {txt}name:  {res}<unnamed>
       {txt}log:  {res}C:\Users\Carolina\Downloads\reuse the seed\Tables 1-4 log_final.smcl
  {txt}log type:  {res}smcl
 {txt}closed on:  {res}29 Jun 2021, 12:04:40
{txt}{.-}
{smcl}
{txt}{sf}{ul off}